


Institutional Archive of the Naval Postgraduate School 


Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations l. Thesis and Dissertation Collection, all items 


1975 


Approximate analytical evaluation of 
extended-Kalman filters. 


Tasdelen, Hasan FOner 


Monterey, California. Naval Postgraduate School 


http://hdl.handle.net/10945/20764 


Downloaded from NPS Archive: Calhoun 


| Calhoun is the Naval Postgraduate School's public access digital repository for 
D U DLEY research materials and institutional publications created by the NP3 community. 
3 Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 
KNOX appointed — and published — scholarly author. 


LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 


APPROXIMATE ANALYTICAL EVALUATION OF 
EXTENDED-KALMAN FILTERS 


Hasan Oner Tasdelen 


DUDLEY KNOX LIBRARY 


NAVAL POSTGRADUATE SCHOOL 
MONTEREY. CALIFORNIA 93940 








NAVAL POSTGRADUATE SCHOO 


Monterey, balifornia 








APPROXIMATE ANALYTICAL EVALUATION 
OF EXTENDED-KALMAN FILTERS 







by 


Hasan Oner Tasdelen 







December 1975 







МТ ез1 Advisor: 


Approved for public release; distribution unlimited. 


1170472 








SECURITY CLASSIFICATION OF THIS PAGE (When Dete Entered) 


REPORT DOCUMENTATION PAGE 


. REPORT NUMBER 2. GOVT ACCESSION NO. 


4. TITLE (and Subtitle) 


Approximate Analytical Evaluation 
of Extended-Kalman Filters 





READ INSTRUCTIONS | 
BEFORE COMPLET:NG FORM 


3. RECIPIENT'S CATALOG NUMBER 















$. TYPE OF REPORT & PERIOO COVERED 
Master's Thesis; 
December 1975 


6. PERFORMING ORG. REPORT NUMBER | 


8. CONTRACT OR GRANT NUMBER(8) 






AUTHOR(®) 






Hasan Oner Tasdelen 













10. PROGRAM ELEMENT, PROJECT, TASK 
AREA & WORK UNIT NUMBERS 





. PERFORMING ORGANIZATION NAME AND ADDRESS 


Naval Postgraduate School 
Monterey, California 93940 






11. CONTROLLING OFFICE NAME AND ADDRESS 12. REPORT DATE 
Naval Postgraduate School December 1975 
Monterey, California 93940 ee 
4. MONITORING AGENCY NAME & AODRESS(I( di(lerent (rom Controlling Office) | 15. SECURITY CLASS. (oí thia report) 
Naval Postgraduate School | 
Monterey, California 93940 UWaecTlas ste 
15a. DECLASSIFICATION/ DOWNGRADING 
SCHEDULE 


16. DISTRIBUTION STATEMENT (of thie Report) 


**Approved for public release; distribution unlimited. 


17. DISTRIBUTION STATEMENT (ol the abetract entered in Block 20, If dillesent from Report) 


- SUPPLEMENTARY NOTES 





19. KEY WORDS (Continue on reveree sida Il necessary and identify by block number) 
Analytic equations, extended-Kalman filter, performance 
evaluation, Monte-Carlo simulation, conditional expectation. 


20. ABSTRACT (Continue on reverse side Il necsssery end identify by block number) 


Analytical equations derived for evaluating linear estimators 
are applied to extended-Kalman filters for approximate perform- 
ance evaluation. Two cases were considered, a single known 
target trajectory and multiple target trajectories with given 
ОТОО NE Sob oeeubrrenee. For the multiple-trajectory case, 
equations are derived for the mean and covariance of estimation 
error in terms of the conditional expectations. Two examples 





FORM 
DD , n 73 1473 EDITION OF 1 NOV 68 iS OBSOLETE 
Pag S/N 0102-014-6601 | Bar 
( age 1) SECURITY CLASSIFICATION OF THIS PAGE (Yhen Dats EntsreJ) 





— o E o s HEP i UU E ERE 
SECURITY CLASSIFICATION OF THIS PAGE(W^hen Deta Entere d. 


are presented to compare the use of the analytical equations 
with Monte-Carlo simulation. 


DD Form 1473 
ean ic 
S/N 0102-014-6601 


SECURITY CLASSIFICATION OF THIS PAGE(When Data Entered) 





Approximate Analytical Evaluation 
ОГ 
Extended-Kalman Filters 


by 


Hasan Öner Tasdelen 
Lieutenant, Turkish Navy 


Submitted in partial fulfillment of the 
requirements for the degree of 


MASTER OF SCIENCE IN ELECTRICAL ENGINEERING 


from the 


NAVAL POSTGRADUATE SCHOOL 
December 1975 





=| í 


En | ; 

n ME E | Ñ 

б? | | 

o - 

š - 
p 
p | - 

| 7 | 
" 
i Ji 
i | | я 


ABSTRACT 


Analytical equations derived for evaluating linear 
estimators are applied to extended-Kalman filters for 
approximate performance evaluation. Two cases were con- 
sidered, a single known target trajectory and multiple 
target trajectories with given probabilities of occurrence. 
For the multiple-trajectory case, equations are derived for 
the mean and covariance of estimation error in terms of 
the сало. Two examples are presented 
to compare the use of the analytical equations with Monte- 


Carlo simulation. 





ША БИЕ ОЕ Б 


fee INTRODUCTION -------------------------------.--- 
II. ANALYTICAL EQUATIONS -------------------------- 
A. ASSUMPTIONS ------------------------------- 

B. DERIVATION OF ANALYTICAL EQUATIONS -------- 


III. APPLICATION OF THE ANALYTICAL EQUATIONS TO THE 
CASE OF MORE THAN ONE TRACK ------------------- 


EEUU US OPBOTATION ------------------- 
B. A LINEAR PROBLEM -------------------------- 
NEM "UATORS FOR NONBINEAR MODELS --------------- 
КЕКТІ: 2-22 --------------- 
Be ТІТКБЕАВТТАТТОМ ----------------------------- 
ОС we Dita FILTER ---------------- 


V. APPLICATION OF THE ANALYTICAL EQUATIONS TO 
EXTENDED-KALMAN FILTERS: A RE-ENTRY PROBLEM --- 


A. THE RE-ENTRY PROBLEM ---------------------- 


B. APPLICATION OF THE ANALYTICAL APPROACH TO 
THE MULTIPLE TRACK NONLINEAR CASE --------- 


C. CONDITIONAL COVARIANCE -------------------- 

VI. A SUBMARINE TRACKING PROBLEM ------------------ 
VII. CONCLUSIONS ----------------------------------- 
APPENDIX A: Computer Programs -------------------- 
A. Program "Eval" ---------------------------- 

B. Monte-Carlo Simulation -------------------- 
COMPUTER PROGRAM ---------------------------------- 
BIBLIOGRAPHY -------------------------------------- 
INITIAL DISTRIBUTION LIST ------------------------ 


10 
12 
12 
14 


26 
27 
22 
50 
LO 
41 
43 


46 
46 





p.l. 
Bac. 


3.3. 
EL. 


25. 


Во. 


Del. 
SO 
5. 
5.4, 
E 5. 
po. 


27: 


EO. 


5.9. 


p.10. 


SINE. 


En 12: 


12. 


BISITEORZRIECURES 


mime MSto eo Mmea posi tilon -= =s- 
БАСЕ не опе velocity ==-============- 
Position estimation error history ----------- 
Velocity estimation error history ----------- 


Variance of position estimation error 


ee з _—_————-------------------- 


Variance of velocity estimation error 
VO. EMO O —— — 


Geometry of Re-entry problem ---------------- 
Altitude history of re-entry vehicle -------- 
Velocity history of re-entry vehicle -------- 
Time history of ballistic coefficient ------- 
Time history of altitude estimation error --- 


Time history of the velocity estimation 
S n СЕС 


Time history of the mean ballistic 
parameter estimation error ------------------ 


Maersnmeevof altitude estimation error vs. 
time ---------------------------------------- 


Variance of velocity estimation error vs. 
BE ME uu... nn 


Variance of ballistic parameter estimation 
Fuss h ss 2. 0. = = — -_-----2------- 


Time history of the mean altitude estimation 
error for 50 Monte-Carlo runs --------------- 


Time history of the mean velocity estimation 
error for 50 Monte-Carlo runs --------------- 


Time history of the mean ballistic parameter 
coto terror for 50 Monte-Carlo runs ---- 


34 


95 
Do 
27 


38 


p 
47 
54 
55 
56 
E 


61 


62 


66 


67 


68 





5.14. 


p.l5. 


E16. 


ERU. 
19. 
EDU. 
20. 


SONA 


Dance of ihe altitude estimation error 


vs. time for 50 Monte-Carlo runs ------------- 


Variance of the velocity estimation error 


VES met su О овор ао Е 


Variance of the ballistic parameter 
estimation error vs. time for 50 Monte-Carlo 


runs ---------------------------------- -- --- 


Mime history of altitude estimation error 


for X(0) = b -------------------------------- 


Time history of mean velocity estimation 


error for X(0) = b -------------------------- 


Time history of the mean ballistic parameter 
estimation error for X(0) =b -------------- 


Eu umce oft altitude estimation error vs. 
time for X(0) = b -------------------------- 


Variance of velocity estimation error 
vs. time for X(0) = b ---------------------- 


5.22. Variance of ballistic parameter estimation 


2. 


5.25. 


о. 25, 


Б 26. 


b 27. 


E25. 


O. 
Өн 25 


error vs. time for X(0) = b ---------------- 


Time history of the mean altitude estimation 
error for two tracks ------------------------ 


Time history of the mean velocity estimation 
error for two tracks ------------------------ 


Time history of the ballistic parameter 
БП тосо ок Шо tracks ============= 


Time history of the variance of the altitude 
estimation error for two tracks ------------- 


Time history of the variance of the velocity 
estimation error for two tracks ------------- 


Time history of the variance of the ballistic 
parameter estimation error for two tracks --- 


R - X filter geometry ----------------- 
cpa cpa 


Tameshistory of the mean Ropa estimation 


27 


78 


86 


87 


88 


89 


90 


91 
93 





ENS. 


6h. 
on 5. 


0.6, 


on 7 - 


EROS 


6.9. 


GNO. 


SEL. 


Time history 


Time history 


Time history 


Variance 


ус. time 


Variance 
vs. time 
Variance 


vs. time 
Variance 
vs. time 


Variance 


mean 6. estimation 


mean 2 estimation error- 10% 





ACKNOWLEDGEMENT 


The author is deeply indebted to the United States 
Navy and the Turkish Navy for providing the opportunity 
to pursue this work. 

I wish to thank my thesis advisor Professor D. E. Kirk 


and all my professors, for their guidance and counsel. 





1. INTRODUCTION 


In state estimation problems one of the tasks is the 
Simulation of a filter or estimator after designing it. 

The purpose of the simulation is to determine the perform- 
ance of the filter under actual operating conditions and 
to investigate sensitivity to inaccuracies or approxima- 
Bons present in the design assumptions. 

Monte-Carlo simulation is a common simulation algorithm 
which is currently used. An important drawback of the 
Monte-Carlo approach is the large amounts of computer 
time generally required to achieve a reasonable degree of 
accuracy. To reduce the computation time requirements, 
fewer Monte-Carlo runs maybe used with an attendant loss 
of accuracy. Because of these limitations of the Monte- 
Carlo approach, an alternative method, cosisting of a 
set of analytical equations, has been derived and a 
computer algorithm has been established for the evaluation 
nma [шло апа prediction error statistics for linear 
filters. It is possible to characterize the propagation 
of the means and covariances of estimation error of a 
filter by difference equations. These difference equations 
can easily be solved using relatively small amounts of 
computer time. The only disadvantage of the analytical 
equation approach is that it only applies to linear filters 


with precomputed gain schedules. 
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ia practice, one may be confronted by nonlinear 
problems, such as space vehicle re-entry or orbit deter- 
mination problems, or fire-control estimation problems. 
A common approach for nonlinear problems is to use an 
extended Kalman filter and to evaluate its performance 
by Monte-Carlo simulation. 

The main objective of this research is to investigate 
the possibility of approximate evaluation of the perform- 
ance of extended Kalman filters by applying the analytical 
equation approach. Іп Chapter II of this thesis the ana- 
lytical equations are derived for linear time-invariant 
estimators. In Chapter III, the application of the ana- 
lytical equations to problems with multiple tracks occurring 
with given probabilities is investigated and an example is 
presented. In Chapter IV, a common filtering approach for 
nonlinear systems is summarized. In Chapter V, the appli- 
cation of the analytical equations to extended Kalman filters 
15 discussed for both the single and multiple track cases 
by using a nonlinear example. Chapter VI contains simulation 
results from another example. The simulations were performed 
using both the analytical equations and the Monte-Carlo 
algorithm for comparison. 

The computer programs that were used are given in 


Appendix A. 
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II. ANALYTICAL EQUATIONS 


Erz ASSUMPT IONS 
In the development of analytes equations the following 
assumptions have been made [1] : 
Ie true State trajectory X(k) for k=0,1,2,.... 
is known. 
2. The measurement equation is 
2(к)=һ(Х(к)) + у(к) (221) 
where 
Z(k) is the q-dimensional measurement vector at 
time t-kT. 
h ШОО unco (which may be nonlinear) of 
X(k). 
V(k) is the measurement noise vector with the 
assumptions 
RS E [ v(x) | UE CIE ГТ 
(The expected value of the measurement 
noise is zero for all time) 
m R(k) for k-j | 
(b). Е[у(юу() | = for all k,j >0 
Q LOr kf j 
(the measurement noise is uncorrelated). 
(c). Since the target trajectory is known 
E[xcvao7] =x(3) .E[yix)] =0 ror all x, jz0 
3, The estimator is linear and described by the equation 


$ cA 8G /-1)*809 [ 2090-68-12] (2.2) 
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mre prediction equation is 
Х(к/к-1)-6 &(k-1/k-1)+A\ U(k-1) (2.3) 
where 
X (k/k) is the n-dimensional estimate of X(k) 
. given measurements 

 0),2(1),...,2(к). 

Х(к/к-1) is the predicted value of X(k) given 
measurements 2(0), 2(1),...,2(к-1). 

g and C are nxn and qxn matrices, respectively, 
and they are assumed to be known. 

U(k-1) is the m-dimensional deterministic forcing 
Meetor at time t-(k-1)T. 

A is an nxm known matrix. 

G(k) is the nxm gain matrix. 

The gains may be found by any method. For example, 
optimal estimation gains which minimize the sum of the 
variance of estimation error can be obtained by using the 


Kalman equations: 


G(k)=B(k/k-1) С? | с мк/к-1) Озң(ю) |74 (2.4) 
р(к/к)= [1-6(к) с | Р(к/к-1) (2.5) 
р(к/к-1)= И Р(к-1/к-1) g7 + Q(k -2) (2.6) 


ША спека лат conditions 

P(0/-1)=Bo =E{[ X(0)-Xo] -[ xco-xoT 

N — 

X(0/-1)=X, 

(In this research gains determined by using the 


Kalman equations have been used for simulations.) 
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Where 
в(к/ю-в4 | к/к) - X00] «[SG/^0 - x09]? | (2.0 
which is the covariance matrix of theoretical ° 


a ci mation error, d 


B(/k-1)-5( [ £G/k-1)-x00] -[£oA-0-x00]7) (2.6) 


Q(k) is the covariance matrix of the random forcing 


ИЕ (К). i.e. 


am Q(k) for k=3 
E [200 197). | (2.9) 
0 Hor К 
If the estimation equation (2.2) is initialized with 
the value 
E1)2X -E x(0) (2.10) 
it can be shown that the optimal estimate (к/к) je Unbiased 
B. e, 
Е| Жж/к) - X(k)] =o (2.11) 
With the assumptions given above one can derive 
difference equations for the mean of estimation error, the 
covariance of estimation error, the mean of N-step prediction 
error and the covariance of N-step prediction error. These 


equations are derived for the time-invariant estimator given 


b euuatrons (2.2) and (2.3) in the following section. 


B. DERIVATION OF ANALYTICAL EQUATIONS 
1. Mean and covariance of estimation error 
The true state of the target at time k is X(k); thus, 
Meta Bauations (2.2) and (2.3) the estimation error X(k) at 


time k is 
14 





><» 


X(k) $ X(k/k) - X(x) 
N(k-1/k-1) *AU(Kel) + ск) 20) 
C 


f X(k-1/k-1) - CAU (k-1)] -X(x) (2.12) 


¿A 


Substituting Equation (2.1) into (2.12) gives 


X(x) =[1 - €(%) с | 8 X(x-1/%-1) + G(%) v(%) 
+ G(x) h(X(x)) - X(0) 


+ [I - G(x) g]A u(x-1) (2.13) 
By defining the deterministic matrices S(k) and D(k) as 
SQ) =[Z - gx) ¢] g (2.14) 


Dix) =[ 1 - GG) cla U(k-1) + Gk) h (х(к)) 
- X(k) | (2.15) 


Equation (2.13) becomes 
K(k) = S(k) X(k-1/k-1) * D(k) + G(x) V(x) (2.16) 
The mean of estimation error is defined as 
m ^ 
Е[ (к)] = Е[ 8(х) #(к-1/к-1) + 200 + G(x) v(x)) 
The matrices S(k) and D(k) are deterministic, hence 
E | a(k) 100) ] = 6(x) E[ Y(0)] = о 
and using properties of the expectation operator gives 
A ^ 
E[X(x)] = so) El £-1/0-1)]) + 200 


Defining 


E[ X(x)] 


Дк) 


gives 


Д(к) = 5(к) Е[ #(к-1/к-1)] + шк) (2.17) 
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mavetion (2.16) can also be written as 


0) Sk) X(k-1/k-1) + D(k) + G(k) V(k) 
+ | S(x) X(k-1) - S(x) X(x-1) | 
= sao X(k-1/x-1) - X(k-1)] AS) 
+ G(k) V(k) + S(k) X(x-1) 
or 
X(k) = S(k) X(k-1) + D(k) + G(k) V(x) 


* S(k) X(k-1) (2.18) 
Since X(k-1) is deterministic and 
E [ Xa-2] ӨЛ. 
Egustion (2.17) can be written as 
А (®) = 8(к)/Д(к-1) + р(к) + 5(к) Х(к-1) (2.19) 
Equation (2.19) defines the mean of estimation error at time 
k, in terms of the mean of estimation error at time (k-1). 


The covariance of estimation error is defined 


as 
P(k/k) = | E) - 400] -| X00) - 400) й 
- E [xao xao? | c 
where 
Y(x) = (к) - (к) (2.21) 
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ne can obtain a 


difference equation for Y(k) by substi- 


ЕПС (2.18) and (2.19) into (2.21), that is, 


1405) < 


Then 


P(k/k) 


Using properties 


P(k/k) 


S(k) [ X(x-1) 


S(k) X(k-1) + 


S(k) X(k- IAS) A Sd 


EA (k-1) - Dk) = S(k 


Z(k-1) 


X(k-1) 
- а-1)| ж а(ю) v(k) 
G(k) V(k) (25 22) 


id 


lI 
(23 


) X(k-1) Y(-1)" 8097] 
) X(k-1) YO" 809^] 
) Y() xG-1) $097] 
O ge)" ] 
of the expectation operator yields 
= S(k) E[ y(k-1) Y(k-1)'| 500)” 
x)" 


+ G(k) R(k) G(k)Í (2.23) 


To obtain expressions for E [ X(k-1) va) ] and E| V (k) Y(k-1)] 


pescan start by 


deriving an expression for Y(0). The 


derinition of Y(0) is 


Ү(0) - 


Х(о) -/Жо) 
Х(0/0) - хо) - Е|Жо)| 
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- X(0/o) - X(0) - E[ £to/0) - xto)] 
Since — E[x(0)] - xto), 
x(0) = X(0/o) - E[ £(0/0)] (2.24) 
But 


X(0/0) = Z(0/-1) + в(о)[ 2(0) - с #(о/-1) ] 


(1- «о) с| “(о/-1) + (о) nato») 
* G(0) V(0) (2.25) 
Substituting (2.25) into (2.24) gives 
x(0) = [1 - eto) e] Жо/-1) + glo) n(xCo)) 
EO) E [L -g(0) с] X(0/-1) 
+ G(O) h(X(0)) + G(0) у 
However, 
E [cto x0] » eco g [vto] 2 o 
е [ с(о) һ0х(0)) ] - e(o) nico) 
and using properties of the expectation operator gives 


убо) = [1 - evo c] $o/-0 * aco atto») 


* g(0) Y(to) - [x - ato) c] E| £(0/-1)] 


- g(o) Artk(o)) 


But Х(о/-1) tarde terministic, known quantity, so 
A 
E[ 2(0/-1)] = X(0/-1) and 


Y(0) = G(0) V(0) (2.26) 
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Nssne Equation (2.22) 
Men (1) Y(0) + G(1) V(1) 
E P C(0) V(0) + G(1) V(1) c 


and it can be shown that 


k-lek-j-1 
үю = ур | G(j) YCJ) +00к) Ук) (2.28) 
o 


Equation (2.28) shows that Y(k) is a linear combination 


of V(k) with CO í clients which are known constant matrices, 


1.е. 
ш =), 1, (к) У(6) (2.29) 
m 
where 
k- 6 -1 
Ly (k) = > Se 602) for 1=0,1,2,..., (21) 
1=0 
(2.30) 


Lk) = G(k) 


From Equation (2.29) it is easily seen that 


Е | Ү(к-1) v(x)? | = E[ v(x) ¥(k-1)7] = 0 
Because of the assumption that E | v(x) v3)" ] = 0 Гог КУ]. 
Ns result reduces Equation (2.23) to 
Р(к/к) = S(k) P(k-1/k-1) $07 
+ G(k) R(k) G(x)? v3) 
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Equation (2.31) is the error covariance propogation 
equation which expresses the covariance of estimation 
error at time k in terms of the covariance of estimation 
error at time (k-1). 
Initial conditions 
To use Equations (2.19) and (2.31) one needs to know 
Mitlal condition values. 


Mie mean of estimation error at time k=O is 


2 (0) 


t 


E [$(o/o) - xo)] 


E [ %(о/0) | - хо) 
Using Equation (2.25) this becomes 


Д4 (оу) = Ef [1 - gto) С) X(0/-1) + G(0) n(X(0)) 


+G(0) V(0)) - X(0) 
But 
E| %(0/-1 ) | = Rey 
Ef c(o) n(x(o))] = &(0) h(xCo)) 
and е | с(о) у(0)| = 0 
thus, 
00) = [ 1600) с] (00/-1) - Хо) 
* G(0) h(x(0)) ( 


Equation (2.32) gives the initial condition for Equation (2.19). 
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The covariance of estimation error at time k=0 is 
Plo/0) = E| X(0) X(0)] (2.33) 
Apesta tuting Equation (2.26) into (2.33) gives 


P(0/0) 


E | с(о) у(о) v(o)l со)” | 


! 


(0) R(o) g(o)" (2.35) 
ЕИ оп (2.35) Бічеѕ ће initial condition for equation (2.31). 


3. Mean and covariance of prediction error 


The one-step prediction is given by 


K(kt1/k) = (к/к) +AU(x) (2.36) 
and the N-step prediction, based on the estimate X(k/k) is 
N-1 
X (een/x) = g x(k/k) + I Oo EA Ui) (2.37) 
1-0 


Defining the deterministic matrix A(N) as 
N-1 | 
AN) =) gT uleti) (2.38) 
і-0 


then 

A ^ 

X(t) = 0 (к/к) + А00) (2.39) 
The N-step prediction error is 


X (krN/k) = X (x+N/k) - X(k+N) (2.40) 


1 


d" (к/к) + A(N) - X(k*N) 


ZI 


E 





The mean of the N-step prediction error is defined as 


S - El 5 (кен/к)) 
= El Y (к/к) + A(N) - X0cN)] 
= ў Е [ £A] + A(N) - X(k#+N) (ОШ, 
But 
Айк) = E| £9. - xao ] 
-Е|Ж(ҡ/к) - x(x) 
and 


Е| Х(к/к)| - Дао + x(x) (2,12) 
Substituting (2.42) into (2.41) gives 

Д.О) = Я Дк) + И (к) + BON) - XOsN) (2.43) 
The covariance of N-step prediction error is 


P(ktN/k) = ЕЕ - (клу) | | блик) 


- кетик) |" (2.44) 
Using (2.40) and (2.43) yields 
ғ en N ^ 
X (k+N/K) Z (k+N/k) = p X(k/k) + T - Х(ұғм) 
2 ИЧ (к) < g X(k) - 


Е [| Х(к/к) = US) 740] 


-g'[Xx^) Am] (2.45) 


+ 


Ze 





which, when substituted into (2.44) gives 


l 


Poen/x) - EV g' [£0 200 ] [Go 


Фо} (и) 
а/ж) (|7 (2.16) 


Equatiows (2.43) and (2.46) give the mean and 


l 


covariance of N-step prediction error based on the mean 
ШОО б Уагтапсе о estimation error at time К. 


In practice, the measurement equations are often 


2(к) = © X(k) + V(x) (2.47) 
Then, in the derivation of the analytical equations, h(X(k)) 
can be replaced with C X(k). This will change Equation 
(2.15) to 
Dik) -[1- $0) || Au(x-1) - хк) (2.48) 
and Equation (2.32) to 
Ao - [1 - (0) 2] [ £(0/-2) - xo] (2.49) 
4, Summary of key equations 
(1). Mean of estimation error 
WE) = SR) (к-1) + (к) + 5(к) X(Oc-1) (2.19) 
(2). Covariance of estimation error 
(к/к) = S(x) Blk-1/k-1) S(k)* 
+ G(k) ROO GR)" (2.31) 
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(3). 


where 


S(k) = [Z - g) c] g (2.14) 
Dik) = [Z - g(x) E ]au(x-1) 
+ GCK)hLXCK)) - X(%) (2.15) 


with initial conditions 
Mio=[1 - ato) e] £to/-1) - xto) 
+G(0) h(X(0)) (2.32) 
P(0/0) = G(0) R(0) а(о)! (9029) 
Mean and covariance of N-step prediction error 


(к/к) = (к) + И хк) + A(N) 


ue» (2.83) 
P(x+N/k) - g" Р(к/к) (Ж x (2.46) 
where 
N-1 | 
БАҚ) - ). gI A ulti) (2.38) 
1-0 


If the system is linear and time-varying the 9,4 


and C matrices are time dependent and are denoted by g(k), 


Z (X), C(k). 


Since 


X(k/k-1) - f(k-1) X(k-1/k-1) +A (k-1) U(k-1) (2.50) 


Z(k) = h(X(k)) + V(k) (2.51) 
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it is only necessary to replace the Z, and Ç matrices with 
@(k-1), Z2 (k-1) and C(k) in Equations (2.14), (2.43), (2.15), 
were) and (2.38). 

The analytical equations derived in this chapter 
are applicable to linear filters with precomputed gain 
schedules. 

From Equation (2.31) it is seen that the covariance 
of estimation error is independent of the track, however, the 


mean of estimation error is track dependent. 
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ШЕШІ APPLICATION OF THE ANALYTICAL EQUATIONS 
TO THE CASE OF MORE THAN ONE TRACK 


The equations derived in the previous chapter are based 
on one.known track. Application of the analytical equations 
using one known track has been studied in Reference / 1 / 
and results have been tabulated. 

IEEE is desired to evaluate the performance of a filter 
for various tracks, a different track can be used for each 
iteration (ensemble member) in a Monte-Carlo simulation. An 
alternative approach based on the previously derived analytical 
equations is: 

(1). Apply the analytical equations to calculate the 
mean and covariance of estimation error for each track indi- 
vidually; these are the conditional means and covariances. 

(2). Calculate the overall mean and covariance using 
the results of (1) and relationships involving the means of 
conditional expectation. 


Assume that there are several tracks to be considered; 
(1) (2) (1) (n) 
О Х(ЮО),.., Х(к) 


where each track has a given probability of occurrence 


Py Por P3520 + Pare Pn ‚ 1.е. 


(1) 
»[ X(k) = х(к) | er 


E a 
py + рә + ..+ р; Е er 
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(1) (i) 

The mean kA (k) and the covariance P(k/k) can be calculated 
(i) (i) 

for the i'th track X(k). Since X (k) has been used for the 
EIA EA 

calculation of the 4 (к), Жұ) is the conditional expectation, 


l.e. 
(i) (i) 

^ c 
Ё) = Е ЕРДЕ - x00 | (3.2) 
(i) 
Which is the mean of estimation error given the track X(k). 
Thus, conditional expectations can be computed for each of 
the tracks. 


To calculate the mean of estimation error, one needs a 


relationship between the conditional means and the mean. 


A. CONDITIONAL EXPECTATION 


The conditional expectation is defined as / 2 /: 
+ OO Qo 


E | A/B=b, | du a fala / B=b,) .da, „да... «da, 


- ОО -DO ers) 


Where 


A is a continuous random vector. 
B is a discrete random vector. 
(а / B-b;) is the probability density function of the 


random vector A given B-b;. 
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The expected value of a continuous random vector A is 


оо + HO 
= М а ГА (а) da, e da, e daa (3.4) 


Equation (3.4) can also be written in terms of the joint 


probability density function of continuous random vectors, i.e. 


+. da dbi... m 


(3.5) 


a fypla,b) . da 


From probability theory the joint probability density function 
can be expressed in terms of conditional density function. 


Since the random vector B is discrete (point conditioning), 


Ш pla, b) 
f (a / B=b.) x ES (3.260) 
J 
or 
fap(2,2) = fa(a / B=b;) - P| 3b; ] (3.7) 


ne Equation (3.7) into (3.5) and replacing the 
integrals with a summation sign for the random vector P 


(since B is a discrete random vector) gives 


Ик ]- Y [sss / B=b,) e da... das 


E. (3.8) 
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The term inside the brackets in Equation (3.8) is the 
conditional expectation given by Equation Wehen 
Equation (3.8) takes the final form of 
B| A ]- 2 ЕГА / B=b; ] P| Beb, | (3.9) 
J 
Equation (3.9) can be used to define the mean of estimation 


error in terms of individual means calculated for each track, 


c 
po - El X(x)] 
B (1) (i) 
- » 360 L кок) | P | x(%)=X(x) ] (3.10) 
= | 
п (i) 
" ) SUP (3.11) 
i-i 
where 


(i) 
ВАЕ P| X(k)=X (x) | 
From Equation (2.31) it is seen that the covariance of 


estimation error is independent of the track for linear 


F Tc s 50 11 сап be calculated once for all tracks. 


B. A LINEAR PROBLEM 


An example of the application of the analytical equations 


to a simplified fire control cstimation problem has been 


29 





presented in / 1 /. In this section, the same problem is 
REumdered Tor three tracks with given probabilities of 
occurrence. The simulation has been done using both the 
analytical equations (program 'EVAL') and a Monte-Carlo 
algorithm. 

The problem is a part of a simplified anti-aircraft fire 
control estimation problem in one dimension. The filter is 
a Kalman filter which is derived by solving Equations (2.2) 
through (2.11) subject to the following assumptions. The 
model of target motion is a point mass moving with constant 


velocity. The state equations are 


X(k+1) = g X(k) + WR) (СОЕ) 
where 
a T 
É = (240) 
0 1 
ze 
ко 
Г- (3.14) 
ü T 


Т іс the time between measurements and equals 1 second. 


The measurement equation 15 


|t 
— 
Si 
N 
lI 


€ X(k) * V(X) 


IN 
ез 
> 
< 

ll 


Ë о] X(x) + v(x) (3.15) 
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where 
X4 (X) ls the range of the target. 
X, (k) is the range rate of the target. 
V(k) is the measurement noise, 
w(k) ls a random forcing input. 
It is assumed that 
g [v0] = 0 (3.16) 
El v(x)?] = 625 (m)? 
Ер И, 
E [w(x)] = o (3.18) 


E [w(x)?] = 225 (m/sec“) © 
= Y (3:19) 


The filter initial conditions are 


EU 102 m 
g[xto)] = (3.20) 
-600 m/sec 
And the initial covariance matrix is 
102 о 
P(0/-1) = 919220 
o 107 


ferme this value for P(0/-1) the gain and covariance 


equations were solved to determine the gain schedule G(k), 


k= ПЕ. 2. ...., 


Sil 





In the simulations three tracks (true trajectories of the 


target) were used with initial conditions 


10), 60 х 102 m 
0) 


-600 m/sec (3.23) 
( 2) 55 x о 
Х(0) = 
-600 m/sec (3.24) 
(3) 70 x 109 m 
Х(0) = 
А 3225) 
-600 m/sec i 


The probabilities of occurrence of the tracks are 


(1) 

Ру = P|X(k)=X(k)] = 0.3 (3.26) 
(1) 

E »[xa)-x ao] - 0.3 o 
(1) 

Р. = P[X(x)=X(k)] = 0.4 (3.28) 


Figures 3.3 throush 3.6 show a comparison of the results 
obtained by using the analytical equations and Monte-Carlo 
simulation. Continuous curves represent Monte-Carlo results 
and triangles represent the results obtained using the 
analytical equations. 10,000 Monte-Carlo runs were used 


(3000 runs for the first track, 3000 runs for the second 
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track and 4000 runs for the third track). From Figures 3.3 
through 3.6 it is seen that the two sets of results are very 


close -- especially the covariances of estimation error. 


D 
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FIG. 3.2. Time history of mean velocity. 
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IV. ESTIMATORS FOR NONLINEAR MODELS 


The extended-Kalman filter is a commonly used approach 
for nonlinear estimation problems. In this chapter, the 
application of extended-Kalman filters to nonlinear models 
is summarized. 

Given a continuous and nonlinear estimation problem, 
in order to apply linear discrete filtering theory the 


problem is first linearized and discretized. 


A. DISCRETIZATION 
Consider a nonlinear continuous system of state and 
observation equations given by 


X(t) ( X(t), U(t), t) (4.1) 


— = 


2(t) 


1 


h(X(t)) + V(t) (4.2) 
The discrete-time representation of these equations has the 
form 
X(k+1) = a(X(k), U(k),k) + W(x) (4.3) 
Z(k) 9 C(X(k) ) + Vk) (4.4) 
Equation (4.3) can be obtained from (4.1) by using the 


relationship 
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/ 


X(k+1) = X(k) * (at): £( X(x), U(k),k ) 


wat = (at)? | Be (X(k), U(k),k) 


l+ 


n n PSS G9, u). 309] 


Ic 


*W(X) CE) 
which is the Taylor series expension of equation (4.1). 
W(k) is a random input which may also be used to approximate 


the effect of the truncated higher-order terms. 


B. LINEARIZATION 

иеге іс available a nominal trajectory X(k) and the 
control U(k) is known, one can linearize equations (4.3) and 
(4,4) by expanding them in a taylor series about the nominal 


trajectory X'(k). For the state equations this gives 


X(x+1) » a( X' (X), U(k),Xk) a 





(оо 54463) 


+ W(k) + Higher order terms (H.O.T) (4.6) 


(X'(k),U(k),k 





Defining the matrix g(k) as 


g(k) - = (X'(k) ,U(k) К) (4.7) 
and truncating the higher-order terms, Equation (4.6) 
reduces to 
X(k+1) = g(k) X(k) + a(X'(x),U(k),x) 
- K(k) X'(k) + W(x) (4.8) 
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Since X'(k) is a known quantity, the second and third terms 
of the right hand side of equation (4.8) are known and 
BEEerninistic. Defining the vector U'(k) as 
U'(k) = al X'(k),U(k),k ) - Z(k) X'(k) (4.9) 
Equation (4.8) becomes 
Х(к+1) = g(k) X(k) + U'(k) + W(k) Comp 
Applying the same procedure to Equation (4.4), the 


linearized form of the observation equation is 


Z(k) = H(k) X(k) + C(X*(k) ) - H(x) X'(k) * V(k) (5.11) 
where 
н(к) - 99 
aX | X'(k) (4.12) 


Defining the vector B(k) as 


B(k) 


A Tk) С) 


yields 


Il 


Z(k) - H(k) X(k) * B(k) * V(k) (4,14) 
Equations (4.10) and (4.14) represent a linear time- 
varying model, U'(k) and B(k) represent bias terms resulting 
from the linearization process. The analytical equations can 

be applied to this model by replacing the U(k-1) term in 


MOS) with the U(k-1) term in Equation (2.37) and defining 


new measurement equations given by 
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Z'(k) = Z(k) - B(k) 


xs OE VO) (4,15) 


C. THE EXTENDED-KALMAN FILTER 
Consider a nonlinear discrete system of state and 


observation equations given by 


X(k+1) = £(X(k),k ) + W(k) (4.16) 
ak) = C( X(k) ) + V(k) (WM) 


In these а ions f and C are nonlinear functions of the 
state variables X(k), W(k) is a random forcing input and V(k) 
is measurement noise with the usual assumptions (an uncor- 
related, zero mean random processes). 


E| w(x) uc | Q(k) 6 


kj 


g|voo vc»? ] - &090 &,. 


in order to apply the linear filter equations, Equations 
(4,16) and (4.17) are expanded about the nominal trajectory 
ЕЕ? if it is available. In practice, it is possible to 
have an idea about the target trajectory for some problems, 
like satellite orbit determination problems, in which case 
Go may predict the target trajectory to be very close to the 
true trajectory and thus, be able to define a nominal 
Eragectory for the purpose of linearization. In such cases 


the gain schedule can be computed before estimation, but in 
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other kinds of problems, it may be impossible to obtain a 
nominal trajectory. An alternative is to linearize the 
problem at each time point about the best estimates of the 
states currently available. In this case the gains can only 
be calculated for each sample as the estimates are available, 
This approach is called the extended-Kalman filter. 

Using an estimate to evaluate the linearization Equation 


(4.16) gives 
х(к+1) = g(k) X(k) + W(x) (4.18) 
where 


à f 
(К) = (4.19) 
2 aX 


X (X/k) 


Similarly Equation (4.17) yields 


Z(k) » H(k) X(k) + V(x) (4.20) 
where 
нк) = 9 р. (4.21) 
ó X X(k/k-1) 


g(x) and H(x) are obtained as indicated in the linearization 
process discussed earlier and evaluated using the estimate as 
the nominal track. The filter estimation equation is 

Х(к/к) = X(k/k-1) + G(k) EG - E(X(x/4-1)) | (4.22) 
The prediction equation is 


X(k/k-1) = Ғ( X(k-1/x-1), kel) (280 


ll 





The gains are obtained by using the relationship 
G(k) = р(к/к-1) НОК)" [ нк) р(к/к-1) НОК)" 
+ R(k) | x 


where 


the theoretical covariance equation is 
(к/к) = [р - ск) нк) | р(к/к-1) 
The covariance propagation equation is 


Е 590) - f(k-1) B(k-1/k-1) g(k-1) + Q(k-1) 


“КОШ he initial conditions are 


X(0/-1) = E [x(o)| 
= X 
m) 
P(0/-1) - P, 


E x0 - X,] [xco - x, | " 





= 


(4.24) 


(4.25) 


(4,26) 
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hee APPLICATION OF THE ANALYTICAL EQUATIONS TO 

EXTENDED-KAIMAN FILTERS: A RE-ENTRY PROBLEM 

In this chapter, the approximate evaluation of the 
performance of an extended-Kalman filter using the analyti- 
cal equations is discussed. Since the analytical equations 
are based on linear estimators with pre-computed gain 
schedules, application of these equations to a problem 
with gains which are not pre-computed is an approximation. 

In an extended-Kalman filter the gains are evaluated on- 
line using the estimates as they are produced. This means 
that the gains will vary from run-to-run even if the track 
is the same. This variation is caused by the differences 
in measurement noise which in turn effect the estimates. То 
use the analytical equation approach to evaluate an extended- 
Kalman filter it is assumed that the linearization which 
results in Equations (4.10) and (4.15) is performed using 
the true Шала This leads to a pre-computed gain schedule 
that is used to approximate the extended-Kalman filter gain 


schedule. 


A. THE RE-ENTRY PROBLEM 


In order to compare the results of using the analytical 
equations and Monte-Carlo simulation, a particular problem 
was selected which contains significant nonlinearities in 
both the state and observation equations. This problem 


deals with estimation of the altitude, velocity and constant 
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ballistic coefficient of a vertically falling body. The 
measurements are taken at discrete time intervals of 

] second by a radar that measures range in the presence 
of (discrete) white gaussian noise. The geometry of the 
problem is illustrated in Figure 5.1. It is assumed that 


the body is falling vertically. 





FIG. 5.1. Geometry of Re-entry problem. 
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The following definitions are used: 


X,(t) +: Altitude 

Xt) : Velocity (downward) 

m Mass (constant) 

Ch : Drag coefficient (constant) 


A Reference area for drag evaluation (constant) 
p : Mass density of the atmosphere 

l. : Radar altitude 

M Horizontal distance 


r(t) EEES 


It is assumed that the effect of gravity 1s negligible. 


The equations of motion are 


Eu. 5 2 (5.1) 
git) = - X4) | 
2m 


The air density is approximated by the exponential function 


D Es CUM (5.2) 
where 
= 5 x 1072 ( 
Defining 
X(t) = E ° (5.4) 


which is a constant ballistic parameter, the state equations 


become 
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EU - -X260 


PS 
ГӘ 
en 
ED 
м2 
| 


PS 

emm. 

ct 

— 
| 


= 0 


which are of the form 
X(+) = £( X(t) ) (5.6) 


The measurement r(t) (range) is given by 
r(t) = V M* « (X,(t) - H )* (5.7) 


and is observed at discrete instants of time so that the 
observed sequence is 
H 2 2 
z(k) = Үл uL tC) ( 5.8) 
where v(k) is white gaussian noise with zero mean and 
constant variance 


E | vce)? | E 


In this ease, the output nonlinearity has the form 
2 2 
cCXQ) ) 2 Va? + (X,(k) - H ) (5.9) 
It is assumed that 
H = 0 


M 


100,000 ft 
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and the filter initial conditions are 


X(0/-1) 


ll 


7 
АС) 


3 x 10? ft 
ц 


DC 
Бех D #71 


mie assumed initial covariance matrix is 


10° О 0 
P(0/-1) = 0 4x108 0 (ocn 
0 0 10 


and the variance of measurement noise is 


Ц 


ft/sec (5.10) 


В = 10.х 10 (гъ) (5.12) 


The true initial conditions of the falling body are 


3x10? ubt 
_ ц 
X(0) = 2x10 ft/sec 
1x10 2 ft! 


Euuns Equation (4.5) to Equation (5.5) it сап be shown 


that the discretized state equations are 


(515) 


X(ktl) = al X(k) ) (5.14) 
KX, (tL) = X,(%) - X,Q0 + 3 X, 00% X,(0) Е 
en 2 ы (куе 
X (ct 1) = Хх (К) - X,(k) X4(X) e j > 
. 


PROD + xy, 092 x (10) Б 
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X (Kk) X,(k) 


(5.16) 





X¿(x+1) = X4(X) (5217) 


Linearizing Equations (5.15) - (5.17) it can be shown that 


the g(k) matrix is 


Mik) = o> ` 
aX X'(k) 
іі "^12. “53 
= от А 22 A23 (5.18) 
A A A 


p 32 39 


where X'(k) is the nominal trajectory and the A,;'s are 


= ' e ? - Yx: (k) 
Ary == (172) X, (k) X3(x) e 1 (5.19) 
A = + X ) KK) -Yxi (k) (5.20) 

A Уе (к) (5.21) 
A15 = M Es (k)“ e y 1 

| ' ОЕ -Ух! (к) 
As 70 X? 007 X33) sx. (k) , с E Y 
- 2YXQ)7 O aN (5.22) 

= ү? | = (CR) 

PP = - 2 X (K) X. (k) e 1 


Cy 00 xx) YX (x) 
2 y X500" X4(x) е 


PSU. X ^ Qe 2 YX, (4) (5.23) 
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yx? 2 2 Ck 

A, = 00)? RK) _ eo X5(k)3 e 29103 
+ 2 X3(k)? ху) 02 YX, (9) (5.2%) 
“31 = 0 (5.25) 
E 70 (5.26) 
A33 ^1 (5.27) 


The H(k) matrix is 
дс 


BOD = SE | x) 


X (k) -H 
үм? + 100 - H)? 
iL (5.28) 
In the simulation of the problem, the true track was generated 
by Using the discretized state equations (5.15) through 
EN The linearization is made about the track in the 
application of the analytical equations and U (k) is defined 


as 


U (k) = a( X'(k),k ) - g(x) X' (xk) (5.29) 


оа 





If the true track is generated by solving the state 


equations then 


X'(k+1) = a( X'(k),k ) (5.30) 
and 

DE - X (k*1) - а) Х (к) (5.31) 

DE) - X (k) - Z(k-1) x'(k-1) (5.32) 


If one uses a track other than that generated by solving the 
state equations, then Equation (5.29) must be used for 
U'(k). The matrices £(x), H(k) and the vector U(k-1) can 
EENUSed Tor application of the analytical equations. In the 
Monte-Carlo simulation of the extended-Kalman filter, one 
must replace X (k) with the estimates (к/к) Equacions 
(5.19) through (5.27) and by X(k/k-1) in Equation (5.28) 
( Blk) must be evaluated about (к/к) and H(k) must be 
evaluated about X(k/k-1) m 

ШІРИгез 5,5 through 5.10 illustrate the results of the 
analytical equations and the Monte-Carlo simulation of the 
extended-Kalman filter. The continuous curves represent 
the Monte-Carlo results (1000 runs) and the triangles 
represent results from the analytical equations. The 
figures show that the analytical equations have predicted 
better performance than that predicted by the Monte-Carlo 
minubeation, Actually, this observation is not generally 


true. The difference between results of the two methods 
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depends on the measurement noise, the nature of the nonlinear- 
ities and the filter initial conditions. For the re-entry 
problem the two results are very close. For example, at 
50,000 ft altitude, the mean altitude error difference 
between the two results is 30 feet. 

In order to investigate the noise dependence of the 
difference between two results, the problem has been solved 
by each method for various covariances of the measurement 
noise and the root mean square of the difference has been 
calculated. The root mean square of the difference between 


two results is defined as 


1 l 1/22 
Ад (1) = EN De | £ (к) ғ & (x) | < (5.33) 
k=1 
N 1072 
| АВ, (қ) - жан) |” 
A VAR I) = == VAR К) - VAR_(k 
el N y | —=] —= 2 (5.34) 
k=1 
where 
Ок), VAR, (x) are the mean and variance of estimation 


error calculated by the Monte-Carlo 
Simulation, 

Ux), VAR,(k) are the mean and variance of estimation 
error calculated by the analytical 
equations, 


N is the number of time points, and 
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Ac (T), A NAR, a (I) are the rms values of the 
difference between the two 


results. 


The results obtained have been tabulated in Table 5.1. 
From, Table 5.1 it is seen that the root mean square 
differences between the two results for both altitude and 
velocity estimation error increase as the covariance of 
measurement noise increases. But this is not true for 
X4 meen is the ballistic coefficient. Table 5.1 indicates 
that estimation of the ballistic coefficient is relatively 
insensitive to the measurement noise level. 

Figures 5.11 through 5.16 represent the results of 
50 Monte-Carlo runs (continuous curves) and the analytical 
equations. It is seen that there is significant difference 


between 1000 Monte-Carlo run-results and 50 Monte-Carlo 


run-results. (See Figure 5.5-5.10). 
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Next, the problem was solved for a different target 


track with initial conditions given by 


3.5 x 102? ft, 


X(0) = О 5 107 ft/sec |= db 85155), 


l x10 1/ғ+ 


The results are illustrated in Figures 517 throush 5.22. 

One can see that the analytical equations predict worse 
performance than that predicted by the Monte-Carlo simulation. 
For comparison, the CPU times used for both method are 

given below. 
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B. APPLICATION OF THE ANALYTICAL APPROACH TO THE MULTIPLE 

TRACK NONLINEAR CASE 

The results derived in Chapter III can also be applied 
to nonlinear problems. Equation (3.11) can be used to 
calculate the mean of estimation error in terms of the 
conditional means of estimation error calculated for each 
track. The conditional means are computed approximately by 
using the analytical equations. 

Since the matrices @(k) and H(k) are calculated using 
the true target trajectory in the nonlinear case, one can 
see that the covariance of estimation error will be dependent 
on the tracks (see Equation (2.31). Thus, one needs a 
relationship to compute the covariance of estimation error 
ie terms of the conditional covariances of estimation error, 
Me, 

(ЗІН C2) (1) (n) 
eee E(k /k), P(K/K),..-,P(k/k),.., P(k/k) ) 
| (5.36) 
iu e o | AU 
where P(k/k) is the conditional е of estimation 
i 


error corresponding to the i'th track X(k). The conditional 


covariance of estimation error is defined as 


(i) m (i) 
Р(К/К) = боуаг | x) /x Q0 -XQ0 (5.37) 
The relationship of the covariance to the conditional 


covariances is derived in the following paragraphs. 
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C. CONDITIONAL COVARIANCE 


Given a discrete random vector B, the expected value of 
the dependent continuous random vector A can be expressed 


in terms of the conditional expectations as 


E[A] 52 Е | А/в=Ъ, | в B=b,] wo 
J 
For simplicity Equation (5.38) can be written as 


E| A ]= E 





E[ A/B] Е (5.39) 


where the subscript B denotes the expectation with respect 
ШО the random vector B. 


mire econditional covariance is defined as [3] 


La/B ” 4 | Шо e [avs] | P - s[2/]]" / ы (5.40) 


[авар (5.41) 


which is the (unconditional) covariance matrix. Replacing 
A with 


[4 - ЕРА | [4 - E ( Al] Е in Equation (5.39) gives ** 


Бл = a 





= | [а-в рај] [а-и В (5.42) 











**Equation (5.39) can be written in the general case as 


Elg(a) y = Ef Bl e(a)/a) | where g(a) may be any function. 
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By adding and subtracting E [ 4/3] inside the parantheses, 


Equation (5.42) can be written as 


I | 7 2] : EA Р oh | 


T 


E Фа] а) з}, 


Р 


ACER 





or 
ЕСЕ - E[A/5]) - Е[А/В /3))* 

I t Е[А/з] у). (Ela T ela)” +@[aa] - [ah 

' @ - Ela/a] Y + ar. Қа) 

ela] - sip] 3) |, (5.43) 
Using the properties of the expectation operator in the 
inner expectation, the cross terms can be calculated as 

E | (a - E[a/8] ) * [ars] - [pta E{(a -E[a/3])/ 3} 


Els/a] - Ela)” 
(5.44) 


(because Е[ А | апа E[ 2/2] are not functions of A, they are 


deterministic with respect to A for a given value of B. 
This means En | Ela/a ] = efa] and E, efaz] /B | = ЕЛ} А/В). 


Using the properties of the expectation operator yields 
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Е{ (А - Е[л/з] у/в}= el2/8] - 4 2l2/8 a} 
= ав] - в[алв]) 
= 0 (5.45) 
Thus, the cross terms drop out in Equation (5.43). Since 
(E [4/2] Е (А) is not a function of A but it is a function 
of B, then Equation (5.43) reduces to 
р, = el Е { (% -Е[а/з]у- (А - Е[А/в] утуз} 
+ Gaa] - [aD - [we] ala)” } 150 


or 


Ра = Е 


ef (a - E[s/3]) - (a - B[a/e] 7/2} | 





G[s/s] -к|А]). [алив] – 3(4])" 


+ 
Е В 


5 








The inner expectation in the first term of the right hand 
side defines the conditional covariance (see Equation (5.40), 


hence 
5,75 (a/a j E { (Е[л/з] - к[А])-(Е[А/в] - [А ү 
(5.48) 


Equation (5.48) gives the required relationship between the 


covariance and conditional covariances, i.e. 
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z (i 
P(k/k) = E] Covar | KG 09 -XG9J | 
X(k) 


(i) (i) Vu 
* E | (к) E ( W(k)- (к) ) 


X(k) 
(D (i) 
= Е } (к/к) | rE) px) - U(x) ) 
X(k) | 
(i) | 
CH Oe) - 4 00 07 | 
X(k) 
which can also be expressed as 
= (4) € CD 
(к/к) = У) Paz), + Уу [2 - йо] 
i-1 i=1 
(1) E 
[ao - Eco]? +, (5.19) 


ЕГІ) (1) 
U (k) and P(k/k) are the conditional mean and 
covariance of estimation error for 
1 
the i'th track X(k), 
Ш(к) ала Р(к/к) are the overall mean and covariance 
of estimation error and 
ls the cite of occurrence of 
i 
ples e I(E) 
First, using Equation (3.11) one can compute the mean of 
estimation error Д(К). Then using Equation (5.49) the 


covariance can be calculated. 
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The re-entry problem was simulated for two tracks 


where 


(1) 
ру = Р {Х(к)=Х(к) } 


! 


0.5 (5.50) 


(2) 
Po = P {X(k)=X(k)} 


0.5 co) 
with 


3 х.102 $5 


(1) 1, 
0) 1 2 хло ft/sec (5.52) 
1 х 1072 1/#% 
3.5 x 10° ft. 
(2) Ц | 
ШШ- | 2 x 10 ft/sec (5 


1 x 10” IV Tr 


Eures 5.23 through 5.28 illustrate the results. The 
continuous curves represent the Monte-Carlo simulation and the 
triangles represent the results from the analytical equations. 
From the figures one can see that the 1000-member ensemble 
Monte-Carlo simulation predicts better performance. Actually 
the differences between the two results are small compared 

to the values of the state variables. For example, the 


maximum difference between the two results is 150 feet for 
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mige altitude error at 300,000 feet altitude (which is the 
initial difference; at later times the two results are 
Closer). The CPU times used for these simulations are 


given below. 


CPU Time 
Сс Сатто simulation (1000 runs).......3 min 50 sec 
Monte Carlo simulation (50 runs)......... O min 17 sec 
Е са салабіопс.....................дпіп у вес 
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VI. A SUBMARINE TRACKING PROBLEM 


This chapter discusses an example with linear state 
equations but nonlinear measurement equations. It is a 
problem in which estimates are made of the states of a 
submerged target: heading, velocity, rest frequency, range 
to the target at the closest point of approach (cpa) and 
distance to cpa. One passive sensor (sono bouy) provides 
measurements of target heading and frequency / 4/. The 
states are 


ста! Range to the target at the closest point 


of approach from the sensor. 
a : XK distance to cpa from the sensor (negative 
before cpa, positive after cpa). 
m врева 
К Target heading. 
Fo : ^ Target rest frequency. 


It is assumed that the target has constant velocity and 
heading. Figure 6.1 illustrates the geometry of the 


problem. The state equations are LA 


Ropa t3) ~ Ropa E о У, % m 
Kopa(Ktl) = X pa(k) + Vo Tr (Ур. % , x) 


Sell 
КО ЕЮ У; ) (6.1) 
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ENG. 6.1. R - X filter geometry. 
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QU = 0 (k) + gl % ) 
FQ) = Е (к) * gut yg) 


with 8, through Es the random forcing terms, hence, 


WO) = ¿(Y o M. % E (6.2) 
S 


Ihe quantities % ү , e are random changes in heading, 
S S O 


velocity and rest frequency. They are assumed to be independent 


Zero mean, piecewise-constant rates of change. They have 


variances defined by 


By lA °] 


S 
2 - gl Y, ? 
5% | 0 | 
I = BL Ye ^] 


The values for the standard deviations were taken from 
typical maneuvering parameters for the target 


Ds HOC) 9/7 1000 see = 1.75533 x 1072 rad/sec 


0; 


Ly = 10 kts / 1000 sec = 5.5555 x 1072 yds/sec* (6.3) 
S 

2e = 0.5 Hz/ 1000 sec = 0.5 x 1072 Hz/sec 
о 


With the expressions for the random forcing terms included, 


the state equation become / Y / 
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Ropa(K*1) = Rn .[K) + A X SO 5 m 


1 
+- = ° EM 4 
X. (k+1) X. Coals ms Tr 2 Ж. 


uncis pr ime > % (m 
v (k+1) =V (k) + I, (m 
(6.4) 
O + Yá ` T 
Ба) = EEC) + Ye T 


where A is +1 for counter-clockwise rotation about the sensor 
Andre tor clockwise rotation about the sensor. A is needed 
to give the correct sign for a given geometry. 


The angle measurement equation is 





RUE) 
0(k) = ө (к) - Tan + | * % (k) for M 0 
: K opal) 
(95) 
A.R pa(k) 
= @(k) - Тат l E -1809 4 у (ку шо Kong? 0 
сра 


The frequency observation equation is 


EK) GV 
ux i. + ҮК) (6.6) 


Y, NEE X opa us) 


2 2 
a ak) T X opa? 


> 





where V. (k), V eK) are measurement noises and N is ihe 


0 
velocity of sound in the medium we is assumed to be 1640 yds/ 


sec). The excitation matrix Q(k) can be found as / 4 / 


Q(k) = E{W(k) ‘W(k)" } 


2— Wl 2. i да, 2 
Ж p. ^ -R е КТ ° : e ; ә 
сра là | cpa cpa là. | E B 2 к 20. i 
B NT. Cm NA EA 1. 
2 2 2 2 2 2 
ТЕУ ЗЕ | T^y 4 ^, T7Y лв Dess MS 
E Y, ера бв | 5 Vo cpa O. | 
A _. ! теи 
Ее: | 
Q(k) = SYMMETRIC | nes 0 | 5 
акы - 
ME | 
TR l6 | 
D. A ON 
| 2 2. 
T. Z 
O 
(6.7) 
The first two terms involve state related terms. 
The state transition matrix is 
1 0 0 0 0 
0 aL T 0 0 
DE IO 0 1 о 0 (6.8) 
0 0 0 1 0 
0 0 0 0 1 
x 
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and the linearized measurement matrix is 


H(k) = 


90(k) Ə0(k) 20(к) Ək) d(x) 


one òX opa ӘҮ o 6. OF, 
эко эю) эку эк) DO) 
Sepa ӘХ ера Е 00, oF, (6.9) 
i A + Xon (k) ER, жез, 
R' (k)? | R'(1)2 
| 
и а с ВАН 
і 
1 2 9 ' 3 | ! ! 
LP) EVE) Ha OO REO VEA) ов! (к) 
Po (Xx) "У, R'(k)? ЩЫ R'(x32 
| | 
O- | I aL 0 
FEM UM T 77 C M E 
f (k) X д E) | pn) | | f (k) 
q ' | 0 
FL) Y, R (K) | EE 
| 
(6.10) 


where X'(k) represents the known track if the analytical 

^ ° 
equations are used and the predicted states (X(k/k-1) ) in 
the actual extended-Kalman filter. ғ (к) is the predicted 


frequency measurement in the extended-Kalman filter and the 





true frequency measurement in the analytical equations. 


R'(k) is defined as 


R'(k) ZEN (к)2 + X pa (0^ (6.11) 


and the covariance of measurement noise matrix is 


2 
о, 0 
R= (6.12) 
7-2 
O Of 
7.61543x 1072 0 


0 Ша ORE 


The filter was simulated using the analytical equations and 
the Monte-Carlo simulation with the initial conditions given 
below, Тһе time between measurements was assumed to be 


100 seconds. 


3.15 х 102 yds. 
-2.61 x 10° yds. 
X(0/-1) 28012039 —vds/sec (6713) 
5.30213 Rad. 
500.432 Hz. 
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ne P(0/-1) matrix is 


1.621 5x107 0 0 0 0 
0 1. 057x109 0 0 0 
Р(0/-1)- 0 0 2.853 0 0 (6.14) 
0 0 0 0.596 0 
0 0 0 O 0.684 


and the true әнмен тг таш conditions are 


2992.59 Yds. 
-5295.85 Yds. 
X(0) = 4,504  Yds/sec. (6.15) 
5.06145 Rad. 


500 Le 


These initial conditions have been obtained from a method 
given in / h /.) 

The analytical equations predict the mean of estimation 
error to be oscillatory. Monte-Carlo simulation results have 
indicated that the filter was unstable. The estimation 
error becomes unbounded as the target passes through CPA. 

The reason is that the system has two angle measurement 


equations used depending on the sign of TA (see Equation 


p 
EU) 3. Since the filter has significant estimation error 
up to CPA, the filter “crosses" the CPA at a different 
measurement time than the track does. Thus, at a certain 


instant the filter uses a different angle prediction equation 


ЭЭ 





then the actual measurement equation. This results in 

a large residual in the estimation equation and the filter 
begins to diverge. If the filter has estimates with very 
small error before CPA, the estimator and the target would 
"cross" the CPA in the same sampling interval. This is a 
very difficult task because the geometry of the system 
indicates that any fluctuations in the target heading largely 


affect the R. and X. 


pa ра” 
In order to force the sign of X (5/1) to change in 
the same sampling interval as the track, frequency measure- 
ments were tested to determine whether they were "up doppler" 
or "down doppler", then the sign of the X o 99 05/41) was reversed 


accordingly so that X. and opal k/k-2) have the same sign. 


pa 
The results of this approach are illustrated in Figures 6.2 
through 6.11 for 1000 Monte-Carlo runs shown as the continuous 
curves. The triangles represent the results obtained from 
the analytical equations. 

Comparing the results of the two methods, one can see 
that they are completely different (especially the covariances). 
There are several possible reasons for this, including: 

(1) The extended-Kalman filter uses the estimates and 
predictions for calculation of the (К), H(k) and Q(k) 
matrices, whereas the true target track was used in the 
analytical equations. Thus, it is reasonable that the gain 


schedules will differ largely if the estimates are poor. 
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(2) Since the linearized measurement matrix H(k) has 
been used, the analytical equations do not include the time 
sharing angle measurement equations; as far as the analytical 
equations are concerned there is only one angle measurement 
equation, because the two equations yield the same derivative. 

(3) Finally, the analytical equations do not "see" the 
change that has been made in the filter algorithm to force 
Me filter to cross the CPA in the same measurement interval 


with the target in the Monte-Carlo simulation. 
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VII. CONCLUSIONS 


It has been shown that the analytical equations can be 
used to evaluate the performance of a linear filter with a 
pre-computed gain schedule. The results obtained from the 
analytical equations are exact and are only closely approxi- 
mated by Monte-Carlo simulations with large ensemble sizes 
and correspondingly large CPU times. 

Por the case of the extended-Kalman filter, it is 
possible to obtain analytical results which approximate 
those obtained by Monte-Carlo simulations. How close the 
results are, depends on the performance of the filter. The 
main difference between the two methods is that one uses the 
true track whereas the other uses the estimates and predictions 
as the nominal trajectory to evaluate the matrices g(k), 

ШОКЕ) (and Q(k),R(k) in some cases, too). If the filter 
performs poorly and gives diverging estimates, one can 
expect that the matrices # (к), H(k), Q(k) and R(k) will be 
completely different in the two methods. This will cause 
different results in the two methods. 

The matrix Q(k) effects only the gain schedule. If it is 
a constant matrix or independent of the states, there will be 
no difference between results obtained by the analytical 
equations and Monte-Carlo simulation. The matrix R(k) has 
two kinds of effects. One is similar to that which Q(k) has, 


the other one is that the R(k) matrix represents the covariance 
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of measurement noise, and in Chapter V, it has been shown 
that the "larger" R(k) is the, larger the difference between 
results becomes. From the results given in Chapter V, it 
is observed that the initial filter states effect the results 
too. Depending on how close the initial filter states are 
to the true initial states, the time at which the filter 
converges will change. 

In summary, one can use the analytical equations to 
approximately predict the performance of an extended-Kalman 
filter. The accuracy of prediction depends on the nature of 


eae problem and the filter. 
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APPENDIX A 
COMPUTER PROGRAMS 


A. PROGRAM "EVAL" 

The Program "EVAL" is a computer algorithm which solves 
the analytical equations derived in Chapter II. The algorithm 
has been written for the special case of the analytical 
equations which is developed for the MK. 86 Fire Control 
System as discussed in Reference / 1/ . The estimation 


equation has been assumed as 


#(к/к) = X(x/x-1) + 800 [ 1,00 200) + 4,00) 2(%-1) 
* MA() 2(k-2) - © &(k/k-1)].— (4.2) 


and the measurement equation as 


A = H XK) 463 Gee) 


where H is the true measurement matrix (C is the measurement 
matrix with synthetic measurements. The matrices M (X), 
M.Ck), M4OO are measurement weighting matrices which have 
been used in the MK. 86 Fire Control System. 

In the development of this research these matrices are 


defined as 


M (K) = 1 (oon identi авга х ) (А.3) 

= = . 
Му (К) = Mo(k) = 0, (А.А) 
Ç =H (A.5) 
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Then Equation (A.1) reduces to the simpler form 


Х(к/к) = #(к/к-1) + о) 2 (0-2 к/к-і))| (А.6) 


With the assumptions given above and starting with the 


definition of the estimation error, i.e. 
X(k) = S(k/k) - X(x) (A.7) 


(where X(k) is the true estimation error vector) 


the following analytical results have been derived in / 1 /: 


WE Covariance of estimation error 
> D c т T 
P(x/x) = S4(k) P(k-1/k-1) 5300)" + 850 ROO S50) 


ТІ: ШІ 
+S, (K) В(к-1) $4007 * $500. R2). 840) 


T T 
tA (k) + A(x)” + A (k) + 8508) (4.8) 


2, Mean of estimation error 


Ë (k) = 3,00) M (x-1) + Dik) + S4(k) X(k-1) (A.9) 


iia lizline values 


P(0/0) = 5 (0) R(0) 8,0) (A.10) 
Й (0) = [1 - 90) с] K(0/-1) - [2 - 90) м,(0)н]х00) 
(А.11) 

4, Mean of N-step prediction error 
Ë. (N+k/k) = d" B) + g^ X(k) - X(Ntk) Ca) 


LES 





5. Covariance of N-step prediction error 


- Ек T 

Бнк) + И" а/о (2) (A.13) 
where 

SQ) 7 gQO 1 (o 


S1 (X) 7 G(X) M, (X) 


S,(k) = G(x) M,(x) 


809 = [1 - 800 2] £ 

A, (%) = 8,00) 5,(%-1) R(%-1) 5, (1)* 

8200 = 8300 [ s,(1-1) 5,(%-2) + 5, (%-1) )] R(x-1) 8,001 
BG) S $ - 1) X(0) + S,(%) H X(%-1) * $509 H X(%-1) 


With the assumptions given by Equations (A.3) - (A.5), 
Equations (A.8) through (A.13) reduce to the form of 
equations derived in Chapter II. 


The computer program given in / 1 / has been used after 
making some small changes and adding three subroutines. The 
program listing is given at the end of this Appendix. In 
order to use the program, one must provide subroutine AK 
which computes the linearized state transition matrix £(k), 
subroutine H(k) which calculates the matrix HKK and 
subroutine UK which calculates the vector Uae ji lle 


matrices Q(k) and R(k) are required to be calculated on-line, 
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then subroutines QON and RON must also be provided. It is 
also necessary to give the proper values to certain flags 


for required operation. 


B. MONTE-CARLO SIMULATION 


A listing of the program used for Monte-Carlo simulation 
is also given at the end of this Appendix. The program has 
been set up for use with both linear and extended-Kalman 
filters. One can simulate a linear or extended-Kalman 
filter by selecting the proper flags and supplying the proper 
subroutines. For extended-Kalman filter simulation, one 
must provide subroutine MEAS for simulation of measurements, 
subroutine TRACK for generating the target track (if not 
read in), subroutine HKK for calculation of H(k), subroutine 
AK for calculation of Y(k), subroutine CK for calculation 
of C(X(k/k-1) ), subroutine XPRED for calculation of 
X(k+1/k) = f ( X(k/k) ), and subroutines QON and RON for 


Q(k) and R(k) if they are not read in. 
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